﻿#include<iostream>
#include<stdlib.h>
#include<stdio.h>
#include<sys/time.h>
#include<omp.h>
#include<arm_neon.h>
using namespace std;

const int N =3000;
const int NUM_THREADS = 8;
bool parallel = true;
float A[N][N];
void Gauss(int n) {
        for (int k = 0; k < n ; k++) {
                for (int j = k+1; j < n; j++) {
                        A[k][j] /= A[k][k];
                }
                A[k][k] = 1.0;
                for (int i = k + 1; i < n; i++) {
                        for (int j = k + 1; j < n; j++) {
                                A[i][j] -= A[i][k] * A[k][j];
                        }
                        A[i][k] = 0;
                }
        }

}
void openMP(){
    int i,j,k;
    float tmp;
    #pragma omp parallel if(parallel), num_threads(NUM_THREADS), private(i, j, k, tmp)
    for( k=0;k<N;k++){
        #pragma omp single
        {
            tmp = A[k][k];
            for(j=k+1;j<N;j++){
                A[k][j] = A[k][j] / tmp;
            }
            A[k][k]=1.0;
        }
        #pragma omp for

        for(i=k+1;i<N;i++){
            tmp = A[i][k];
            for( j=k+1; j<N;j++){
                A[i][j] = A[i][j] - tmp * A[k][j];
            }
            A[i][k]=0.0;
        }
    }
}

void openMP_simd(){
    int i,j,k;
    float32x4_t vt;
    #pragma omp parallel if(parallel), num_threads(NUM_THREADS), private(i, j, k,vt)
    for( k=0;k<N;k++){
        #pragma omp single
        {

            vt = vmovq_n_f32(A[k][k]);
            for (int j = k + 1; j < N; j+=4) {
                if (j + 3 >= N) {
                    while (j < N) {
                        A[k][j] /= A[k][k];
                        j++;
                    }
                }
                else{
                    float32x4_t var = vld1q_f32(A[k] + j);
                    var = vmulq_f32(var, vt);
                    vst1q_f32(A[k] + j, var);
                }
                A[k][k] = 1.0;
            }
       }
            #pragma omp for
            for(i=k+1;i<N;i++){
                //tmp = A[i][k];
                for( j=k+1; j+3<N;j+=4){
                    if (j + 3 > N) {
                        while (j <= N) {
                            A[i][j] -= A[i][k] * A[k][j];
                            j++;
                        }
                    }
                    else {
                     float32x4_t v1 = vld1q_f32(A[i] + j);
                     float32x4_t v2 = vld1q_f32(A[k] + j);
                     float32x4_t v3 = vld1q_f32(A[i] + k);
                     v2 = vmulq_f32(v3, v2);
                     v1 = vsubq_f32(v1, v2);
                     vst1q_f32(A[i] + j, v1);
                    }
                }
                A[i][k]=0.0;
            }

    }
}
int main()
{
        double time=0.0;
        for (int i = 0; i < N; i++)
        {
                for (int j = 0; j < i; j++)
                        A[i][j] = 0;
                A[i][i] = 1.0;
                for (int j = i + 1; j < N; j++)
                        A[i][j] = rand();
        }
        for (int k = 0; k < N; k++)
                for (int i = k + 1; i < N; i++)
                        for (int j = 0; j < N; j++)
                                A[i][j] += A[k][j];

        struct  timeval  tv_begin,tv_end;
    gettimeofday(&tv_begin,NULL);
    //openMP();
    //Gauss(N);
    openMP_simd();
    gettimeofday(&tv_end,NULL);
    time+=(tv_end.tv_sec-tv_begin.tv_sec)*1000000+(tv_end.tv_usec-tv_begin.tv_usec);
    cout << time<<endl;

    return 0;
}

